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The q-model is a random walk model used to describe the flow of stress in a stationary granular 
medium. Here we derive the exact horizontal and vertical correlation functions for the q-model in 
two dimensions. We show that close to a critical point identified in earlier work these correlation 
functions have a universal scaling form reminiscent of thermodynamic critical phenomena. We 
determine the form of the universal scaling function and the associated critical exponents u and z. 



I. INTRODUCTION 

The q-model was introduced by Coppersmith et al. pi 
to describe the flow of stress in a static granular medium. 
Although it does not provide a complete solution to the 
problem of stress flow in a granular medium, it provides 
a good flrst approximation. Moreover, the q-model is 
closely related to models that describe the process of ag- 
gregation in statistical mechanics '5], the transport of 
electrons on the surface of an integer quantum Hall mul- 
tilayer [5 , passive scalar turbulence , and the branch- 
ing of river networks [5], among others. In earlier work 
by Lewandowska et al. [Sj it was shown that although the 
q-model describes physics far from equilibrium, nonethe- 
less, for a particular value of its parameters, the q-model 
behaves in a manner reminiscent of the critical point 
in an equilibrium thermodynamic phase transition. It 
was found that the q-model is remarkably soluble and a 
number of exact critical exponents and universal scaling 
functions near to the critical point were obtained. The 
purpose of the present work is to strengthen the analogy 
to thermodynamic critical phenomena by calculating the 
two point correlation functions for the q-model. 

In experiments that the q-modcl is used to describe a 
pack of beads is loaded from above with a uniform stress. 
The distribution of load at the bottom of the bead pack 
as well as the propagation of stress through the pack 
have been measured experimentally 0. In the q-model 
it is assumed that the beads lie at the sites of a regular 
lattice. The beads in each layer are supported by their 
nearest neighbors in the layer below. For simplicity let us 
consider the q-modcl in two spatial dimensions. Let us 
suppose that the beads are arranged in a square lattice 
as shown in the flgure. For each bead a fraction / of its 
load is supported by its neighbor to the left in the layer 
below and a fraction 1 — / by its neighbor to the right. 
These fractions are assumed to vary randomly from bead 
to bead. We denote the depth of the layer by t and the 
position of a bead within a layer by n. The load on 
bead n in layer t -I- 1 is then determined by the recursion 
relation 

Wn{t + 1) = fn+l,tWn+lit) + (1 - /«-l,t)u;„_i,t (1) 

Given the random fractions / one can use eq ([!]) to propa- 
gate an applied uniform load in the top layer downwards. 
The goal is to obtain a statistical description of the load 



as a function of depth by analytic or numerical solution 
of eq 

In the q-model it is generally assumed that the frac- 
tions are independent identically distributed random 
variables with a distribution P{f) that satisfies the sym- 
metry requirement P{f) — P{1 — f). The simplest choice 
is to assume that each fraction is either zero or one with 
equal probability; this is called the singular distribution 
and it was shown in ref [6] that the behavior of the q- 
model for this distribution is reminiscent of a critical 
point in equilibrium thermodynamics. Other distribu- 
tions can be characterized by a parameter 5 (defined pre- 
cisely below) that measures how far the distribution is 
from the critical point. 

In this paper we calculate the horizontal and verti- 
cal two point correlations of the load on the beads by 
generalizing the methods of ref [6]. The horizontal load 
correlation 

C{m,t) ^ {Wn{t)Wn+2m{t)) (2) 

is the covariance of the load on two beads that are in 
the same layer at depth t and have a horizontal sepa- 
ration 2m. The brackets (. . .) denote an average over 
the random fractions. By translational invariance the co- 
variance depends only on the separation of the two beads 
(2m) and not on the absolute location of either bead (n 
and n + 2m respectively) . By analogy to thermodynamic 
critical phenomena we would expect the horizontal cor- 
relation length to diverge as ^horizontal 1/5'^ and the 
vertical correlation length to diverge as ^vertical fx 1/5'^'^ 
close to the critical point. Furthermore we would expect 
the horizontal correlation function to have the universal 
scaling form 

Ch{m,t)-l^^g{m5''Mn- (3) 

We find that in two dimensions the correlation function 
does indeed have a universal scaling form for large depth 
t and small 5 that we determine. Moreover we find the 
correlation length exponent v = 1, the dynamical expo- 
nent z = 2 and the amplitude exponent a = 0. 
The vertical two point correlation 

C„(m, t, t) = {Wn{t)Wn+2m{t + 2t)) (4) 

is the covariance of the load on two beads that are lo- 
cated in horizontal layers separated by a depth of 2r and 
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FIG. 1: The q-model in two dimensions (reproduced from 
ref |6]). The beads are assumed to be arranged in a square 
lattice. Each bead is assumed to be supported by its near- 
est neighbors in the layer below. A random fraction / of 
the load is transferred to the neighbor at left; the remainder, 
(1 — /) to the neighbor to the right. The fractions vary from 
bead to bead and are assumed to be independent identically 
distributed random variables, t specifies the horizontal layer 
occupied by a bead. The horizontal layers are numbered by 
consecutive integers. The position of a bead within a hori- 
zontal layer is given by n. n is even in a layer with even t; 
odd, for odd t. Thus the horizontal separation of two beads 
in the same layer must be an even integer. Similarly the hori- 
zontal separation for two beads in layers with an even vertical 
separation must also be an even integer. 

that have a horizontal separation of 2m. Since the q- 
model is only translationaUy invariant in the horizontal 
direction the vertical two point correlation depends on 
both the separation of the two layers (2t) as well as the 
absolute location of the first layer (t) . Again by analogy 
to thermodynamic critical phenomena we would expect 
the vertical correlation to have the universal scaling form 

c,(TO,i,r) - 1 = ^H(m<5^^<5''^T^''"). (5) 

We find that in two dimensions the correlation function 
does indeed have a universal scaling form for large depth 
t and small S that we determine. Moreover we find the 
correlation length exponent 1^=1 and the dynamical ex- 
ponent z = 2 in agreement with the values obtained from 
the horizontal correlation calculation and the amplitude 
exponent 13 = 1. 



The relationship of the q-model to random walk mod- 
els can be understood as follows. Suppose that instead of 
a uniform load, a unit load is applied to a small number 
of beads in the top layer. In the singular case the sub- 
sequent propagation of the load follows the trajectories 
of random walkers that coalesce upon contact and move 
together thereafter. The non-singular case can be in- 
terpreted as walkers that have a certain probability to 
fission. By contrast in common random walk models 
the walkers either pass through each other in the non- 
interacting case or bounce off each other or annihilate 
each other upon contact in interacting models [85 . It will 
be seen below that the horizontal correlation functions of 
the q-model decay as Gaussians rather than exponentials 
reflecting their close kinship to random walks. 

In the prior work of Lewandowska et al. [5] the critical 
point was characterized by computing the covariance of 
the load on a single bead. Ref [6] derived the univer- 
sal scaling function that describes the evolution of this 
quantity and computed a number of critical exponents 
including the product vz = 2 (but not the separate val- 
ues of v and z). The present work adds to this body of 
knowledge the universal scaling forms for the two point 
correlations and the important critical exponents v and z. 
The present work is limited to two dimensions where the 
critical behavior is most interesting. Ref [6^ also studied 
higher dimensions and thereby discovered that the upper 
critical dimension for the q-model is three. Moreover ref 
considered the effect of an "injection term" that takes 
into account the weight of the beads which has been ne- 
glected in comparison to the applied load in eq ([T]). In 
other related work Rajesh and Majumdar 9J computed 
the two point correlations for the q-model but at large 
depth and far from the critical regime. Finally Snoeijer 
and van Leeuwen refs HUl [H] have studied the distribu- 
tion of load at asymptotically large depths, as discussed 
further in the conclusion. 



II. HORIZONTAL CORRELATIONS 

A. Exact Solution 

Our purpose in this subsection is to evaluate the evolu- 
tion of the correlations eq ([2| with depth by use of eq ([I]) . 
We begin by noting that the random fractions that ap- 
pear in eq ([T]) are assumed to be independent identically 
distributed random variables with a symmetric distribu- 
tion P{f) that satisfies the requirement P{f) = P(l — /). 
The distribution of fractions may be characterized by a 
parameter e defined via 

(6) 

It is easy to verify that e = 1 for the singular distribution 
and e = I for the uniform distribution. Thus 5 defined 
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as 



(5 = 1 -e 



(7) 



is a measure of the distance of a distribution from the 
singular distribution. 

We assume that a uniform load is applied to the top 
layer. Hence the correlation in the top layer Cm(?n, 0) = 
1. As shown in section II A of ref 6J, the subsequent 
evolution of the correlation is governed by 



Ch{m,t + 1) = y^Hmnc{n,t) 



^ - ^^rnS ] Ch{m - l,t) 



(8) 



The evolution matrix H may be interpreted as the Hamil- 
tonian of a quantum particle on a lattice with a non- 
Hermitian barrier at the origin. It was shown in section 
II B of ref ^ that one can obtain an integral expression 
for the correlation at depth t by making a bi-orthogonal 
expansion in the left and right eigenstates of the non- 
hermitian Hamiltonian H . An important subtlety in this 
procedure is that one has to verify that the eigenstates 
of H are complete since there is in general no guarantee 
that a non-Hermitian operator will have a complete set of 
eigenstates. Making use of the bi-orthogonal expansion 
and the expressions for the left and right eigenvectors of 
H derived in ref [3] we obtain 



Ch{m,t) = 1 + 



dk 



cos 



2tt V2, 
[—ie^ sinfc — e(l — e)(l + cosfc)] 
[(l-2e + 2e2) + (l-2e) cos k] 



(9) 



Eq ^ is an exact expression for the horizontal correla- 
tion function and is the main result of this subsection. 
We note that eq (|9| applies for all m except m = 0. 
For TO = the horizontal correlation function Ch(0,t) co- 
incides with the variance calculated by ref [S] and the 
expression for c/j(0,t) is given by their eq (41). 



B. Critical Exponents and Scaling Limit 

To obtain the scaling limit of the correlation function 
we now simplify eq Q, assuming t ^ 1 and 5 1, but 
without making any assumptions about the relative size 
oft S and toUSI. We find 



Ch{m,t) - 1 = -}[0,ii) 



d_ 



/(^,A*)) (10) 



where 6 = is a scaled measure of depth and /z = 
2mS is a scaled measure of horizontal separation and the 
function 



du 



1 



(11) 



It is evident from eqs (10) and (11) that the horizontal 
correlations have the expected scaling form eq ([S]). Fur- 
thermore, from the form of the scaled variables 9 and /i 
we infer that the horizontal correlation length exponents 
V = 1 and the dynamical exponent z = 2. 

To derive the asymptotic behavior of the correlation 
function in the scaling regime it is convenient to rewrite 
the correlation function as 



Cfi{m,t) — 1 = — f dsexp ~^ ("^ ^ ^) ^"^P ^^4'*^ 

(12) 

This form is obtained by noting that f(9, fj,) is an in- 
tegral over a product of a Gaussian and a Lorentzian. 
Using Parseval's theorem it may therefore be written as 
an integral over a product of the Fourier transforms of 
the Gaussian and Lorentzian factors. Using this trans- 
formed representation for f{0,fi) leads from eq (10 1 to 
eq (fT2l). 



Now let us consider the regime of small depth, 6 <C 
1 or equivalently tS"^ <C 1 (the "critical regime"). In 
this regime we would expect the correlation function to 
be indistinguishable from the correlation function at the 
critical point. Indeed we find that for small 9 eq (12) 
simplifies to 



Chim,t) - 1 



1 



ds exp 



1 



(13) 



Thus the correlation function depends only on 11/ 9 or 
equivalently mj^/i, and is independent of 8. Eq (13) re- 



veals that in the critical regime the loads (or to be precise 
the deviations in load from the mean) on neighboring 
beads are s tron gly anti- correlated. For small distances 
(m/y/i) eq (13) simplifies to 



Ch{m,t) - 1 



1 



2to 



(14) 



At large distances (m/Vt ^ 1) eq (13) simplifies to 

,2 ^ 



Ch{m,t) - 1 



t 



exp 



m 

T 



(15) 



in other words the anti-correlations decay as a Gaus- 
sian. In summary we find that in the critical regime 
there are strong horizontal anti-correlations whose range 
grows with the square root of the depth. 

Next let us consider the regime of large depth, 9^1 
or equivalently tS"^ 3> 1 (the "saturated regime"). In this 
regime eq (12) simplifies to 



Ch{m,t) = - 



1 



exp 



TO 

T 



(16) 
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Thus in the saturated regime there are weak anti- 
correlations. The amplitude of these anti-correlations 
falls off inversely with depth whilst their range contin- 
ues to grow as the square root of the depth. 



III. VERTICAL CORRELATIONS 

A. Exact Solution 

There is a remarkably simple exact relationship be- 
tween the horizontal and vertical correlations for the q- 
model [9]. In order to establish this relationship it is in- 
structive to first consider the correlations between beads 
in two consecutive layers, t and t + I, with a horizontal 
separation 2m + l. By use of the evolution eq ([T]) we may 
write 

{Wn{t)Wn+2m.+ l{t + I)) = (/n+2m+2,t W„ (t) W„+2m+2 (<)> 
+ {{1 - fn+2m,t)Wnit)Wn+2mit)) 

(17) 

Note that the load on the beads in layer t is independent 
of the fractions that appear on the right hand side of eq 
(171; these fractions determine the subsequent propaga- 



tion of load from layer t propagates to layer t+1. Hence 
the averages on the right hand side of eq (17) factorize. 



Recalling that by virtue of the symmetry condition on the 
distribution of the fractions, (/) = ((1 — /)) — 1/2, and 
making use of the definition of the horizontal correlation 
eq ([2]) we obtain 

{Wn{t)Wn+2m+l{t + I)) = ^Ch{2m + 2, t) + ^Ch{2m, t) . 

(18) 

Similarly one can relate the correlation between beads 
separated by two layers, Cy{m,t,T — >■ 1), to horizontal 
correlations by twice using the evolution eq ([T]) to obtain 

Cy{m,t, 1) ^ ^Ch{m + 2,t) + ^Ch{m,t) + ^Ch{m - 2,t). 

(19) 

By now the astute reader will have noted the appearance 
of binomial coefficients in eqs (flSl) and ( 19 1 and indeed 



it is not difficult to show by induction that 
, . . 1 / 2t 



Ch{m-k,t). (20) 



Eq (20 1 is the exact relationship between vertical and 



horizontal correlations and is the main result of this sub- 
section. 



B. Critical Exponents and Scaling Limit 

In this section wc derive the universal scaling behavior 
of the vertical correlation function. There are three cir- 
cumstances to consider. It is clear from fig [l] that each 



bead sits at the bottom of an inverted cone of beads 
that are partially supported by it. The correlation be- 
tween two beads will evidently be strong if the upper 
bead of the pair lies inside the support cone of the lower 
bead. Under this circumstance, which wc call the "di- 
rect case", |m| < r. On the other hand if the two beads 
are so far apart that their support cones do not intersect 
at all the correlation between their loads is rigorously 
zero. The "uncorrelated case" corresponds to the con- 
dition \m\ > 2t + T. Finally there is the intermediate 
"indirect case" in which the two support cones do inter- 
sect but the upper bead does not lie inside the support 
cone of the lower bead. 

We first consider the scaling limit for the direct case. 
The sum in eq (20) can be separated into two parts. The 



first part corresponds to the single term with k — m. 
This term determines the contribution of the covariance 
in load of the upper bead to the vertical correlation with 
the lower bead. It will emerge that this co-variance con- 
tribution leads to positive correlations and tends to dom- 
inate the direct vertical correlations. The second part 
corresponds to the remaining terms in eq (20). It will 



emerge that this contribution leads to anti-correlations 
and is generally subdominant. 

We now compute the co-variance contribution to the 
direct vertical correlation. Retaining only the k — m 



term in eq ( 20 1 we obtain 



Cy{m,t,T) = 



22t 



2t 



(21) 



Assuming that r 3> 1 and |m| 3> 1 we may approximate 
the binomial coefficient as a Gaussian. Finally using the 
scaling limit of Ch{Q, t) given in eq (48) of ref |6| we obtain 



Cy{m,t,T) - 1 = 



: exp 



{msy 



1 -j= I ds exp(— sv t(5^)e 



(22) 



Eq ( 22 ) gives the covariance contribution to the the direct 
vertical correlations in the scaling limit. It agrees with 
the conjectured universal scaling form in eq ^ and we 
find that the critical exponents are /3 = 0,z^ = l,2; = 2. 

To gain further insight into the covariance contribution 
to the direct vertical correlations it is useful to consider 



the limiting behavior of eq ( 22 ) in the critical regime 
[tS"^ ^1) and the saturated regime [15"^ ^1). In the 
former limit the dependence on 5 cancels and we find 



Cy{m, t, t) 



exp 



m 

T 



(23) 



Thus in the critical regime there are strong vertical corre- 
lations (not anti-correlations) and these correlations grow 
with the square root of the depth t. The horizontal range 
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of the correlations grows as the square root of the verti- 
cal separation between the two beads, r. If we imagine 
increasing the horizontal separation of the beads, keep- 
ing the vertical separation fixed, then the correlations die 
away well before we reach the boundary of the support 
cone. In the saturated regime we find 



1 1 



cxp 



m 

T 



(24) 



Here too there are strong vertical correlations (not anti- 
correlations) but these correlations have saturated and 
are no longer growing with depth t. Their horizontal 
range grows as the square root of the vertical separation 
between the beads, r. The correlations die away well 
before we reach the boundary of the support cone, just 
as they do in the critical regime. 

We turn now to the second contribution to the verti- 
cal correlations in the direct case, namely the sum over 
terms with fc ^ m in eq ( 20 ) . In the scaling regime we 



may approximate the binomial coefficients in eq ( 20 1 as 



a Gaussian, convert the sum into an integral and extend 
the range of integration to infinity to obtain 



Cy{m,t,T) = 



dk 



1 



cxp 



c/i(m -k,t). 



(25) 

Using eq ( 12 ) for Ch wc may perform the integral over k 



explicitly to obtain 

Cy{m,t,T) - 1 = - 




exp 



Note that this contribution to the vertical correlation 
function is also of the conjectured form eq ([s]) with the 
same critical exponents /3 = O^v = 1 and 2 = 2. It is 
noteworthy that this contribution is a function only of 
the {t + t)6^ and the ratio md / ^/{t^^T)P ; it is a more 
restricted function of its arguments than required by the 
scaling form eq ([s]). The full vertical correlation in the 
direct case is the sum of eq (22 ) and eq (p6|); however, as 



noted above and shown below, the contribution of eq ( 26 ) 
is generally dominated by the contribution of eq (22 1. 



To gain further insight into the second contribution 
to the direct vertical correlations it is useful to consider 



the limiting behavior of eq ( 26 ) in the critical regime 
{t6^ <C 1) and the saturated regime {tS'^ ^1). In the 
former limit the dependence on S cancels and we find 

Cy (m, t,T) — 1 = exp / as e ' 

\ t + Tj Jo 



X exp 



(27) 



Eq (27) is still a formidable expression but it simpli- 
fies to Cy{m,tT) — 1 sa —1 for m/^/t + r ^ 1 and to 



Cy{'m, t) — 1 — 2 for m/\/t + t 1. Thus the second 
contribution to the direct vertical correlation in the crit- 
ical regime corresponds to a weak anti-correlation that 
is negligible compared to the strong positive correlation 
implied by the first contribution, eq (23). Next let us 



consider the saturated regime {t6^ > ly In this regime 
eq ( 26 ) simplifies to 



Cy{m,,t,T) 



1 1 

I — — , exp 



m 

't + T 



(28) 



Eq ( 28 1 should be compared to the first contribution eq 



( 24 ). Thus in the saturated regime too the second contri- 



bution to the direct vertical correlations is smaller than 
the first contribution. 

In summary, in the direct case the vertical correlation 
is the sum of two contributions both of which have the 
conjectured universal scaling form, eq ([s]), with expo- 
nents /? = 0, 1/ = 1 and z — 2. The first contribution eq 
( [22| ) corresponds to a positive correlation and generally 
dominates. The second contribution eq ( |26[ ) corresponds 
to an anti-correlation and is generally subdominant. It is 
illuminating to look at the limiting behavior of the two 
components in the critical regime {tS'^ <C 1) and the satu- 
rated regime. In the critical regime the first contribution 
is given by eq (23) and the second contribution, which 
lies between —1 and —2, is given by eq (27). In the sat- 



urated regime the first contribution is given by eq (24) 
and the second contribution by eq ( 28 ) . 



Finally, we now turn to the vertical correlation in the 
indirect case, |m| > r. In this case the first contribution 
is absent since the term k = m lies outside the range of 



the sum in eq ( 20 ) . Thus in the indirect case there is only 



a weak anti-correlation between beads given by eq ( 26 1 



This concludes our analysis of the vertical correlations. 



IV. CONCLUSION 

In this paper we add to the existing body of known 
results about the q-model by obtaining new results on 
the exact horizontal and vertical correlation functions. 
We show that the correlation functions have a universal 
scaling form close to the critical point identified in pre- 
vious work [S] and we determine the scaling form and 
associated critical exponents. It is intriguing that the 
q-model, which describes physics far from equilibrium, 
nonetheless shows behavior reminiscent of equilibrium 
thermodynamic critical phenomena. Among the many 
questions that our work leaves open we here mention two. 
First we note that it may be possible to calculate three 
point and higher correlation functions by means of a non- 
hermitian generalization of Bethe ansatz [12j. Second, it 
is desirable to determine the dynamics of the entire dis- 
tribution of load on a single bead. Ref obtained the 
dynamics of the entire distribution right at the critical 
point and Snoeijer and van Leewen llj studied the large 
depth asymptotic dynamics of the load distribution in a 
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case far from the the critical point. Close to the critical 
point ref made a scaling hypothesis about the form of 
this distribution and it remains of interest to verify this 
conjecture via numerics or exact solution. 
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First, approximate cos^*(fc/2) « exp(— fc'^t/4). This form 
is justified for small k and leads to negligible error for 
large k since both the exact expression and the approxi- 
mation are negligible for k large compared to Xj^/t. For 
the same reason we may extend the range of integration 
in eq (|9| to infinity with negligible error. Finally since the 
integrand has negligible weight for large k we may expand 
both numerator and denominator in the secon d lin e of eq 
([9| to leading order in k and 5. Eqs ( 10 \ and (111 result 



after making these approximations and a rescaling of the 
integration variable. 

The alert reader will note that we have also transformed 
eq (48) of ref 16! , which features an integral of the product 
of a Gaussian and a Lorentzian, into an integral over 
the corresponding Fourier transforms instead by use of 
Parseval's theorem. 



